	program ProbilityDensity

	implicit none

	integer :: datwid,i,j,k,nargs,ios,steps,column,Ndat
	real :: datmin,datmax,difdat,tmp(26),p(1000)
        !parameter(datwid=100,steps=1500000)
        parameter(steps=2500000)
	character(len=200) :: arg,datfile,outfile
	integer :: hist(1000)
	real :: data(steps)

        nargs = iargc()
        if(nargs /= 4)then
            write(6,*) "the number of argument is incorrect"
            write(6,*) "usage: $0 datfile 2 100 PD.dat"
            goto 1000
        endif

      	call getarg(1,arg)
      	datfile = arg
        call getarg(2,arg)
        read(arg,'(i4)')column
        call getarg(3,arg)
        read(arg,'(i4)')datwid
        call getarg(4,arg)
        outfile = arg

      	open(10,file=datfile,status='old')
      	open(20,file=outfile,status='unknown')
	i=0
        datmin=1.0E10
        datmax=-1.0E10

        readdata: do
            read(10,*,iostat=ios)(tmp(j),j=1,column)
            if(ios .ne. 0) exit readdata
            i=i+1
            data(i)= tmp(column)
            datmin = min(datmin,data(i))
            datmax = max(datmax,data(i))
	end do readdata
	Ndat=i
	write(*,*) 'There is ',Ndat,' lines.'

        do i=0,datwid
            hist(i)= 0
        end do
        difdat  = datmax - datmin
        do j=1,Ndat
            i=int(((data(j)-datmin)/difdat)*datwid)
            hist(i)=hist(i)+1
!	    write(*,'(I5,F5.2,1X,I5,1X,I4)') j,data(j),i,hist(i)
        enddo
        do i=0,datwid
            p(i)=hist(i)*1.000/Ndat
!	    write(*,*) i, p(i)
	    write(20,'(F7.3,1X,F7.4)') (i*1.0/datwid)*difdat+datmin,p(i)
        end do

	close(10)
	close(20)
1000	end program
